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Abstract 


This document presents the mathematical structure of the programs written to construct a 
nonlinear predictive model to forecast solar flux directly from its time series without refer- 
ence to any underlying solar physics. This method and the programs are written so that one 
could apply the s am e technique to forecast other chaotic time series, such as geomagnetic 
data, attitude and orbit data, and even financial indexes and stock market data. 

Perhaps the most important application of this technique to flight dynamics is to model God- 
dard Trajectory Determination System (GTDS) output of residues between observed position 
of spacecraft and calculated positon with no drag (Drag Flag = off). This would result in a 
new model of drag working directly from observed data. 
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Section 1. Introduction 


This document presents numerical techniques for constructing nonlinear predictive models 
to forecast solar flux directly from its time series. As a continuation of our previous 
research in unders tandin g the dynamics of solar activity (References 1 through 5), we 
consider the dynamic evolution of our system (solar activity) in a re-constructed phase 
space that captures the (strange) attractor,* and we give a procedure for constructing 
parameterized maps that describe the evolution of points in the phase space into the 
future. The predictor would necessarily depend on past data points and different itera- 
tions of the map. Thus the map is regarded as a dynamical system and not just a fit to 
the data. The invariants of our dynamical system, the Lyapunov exponents and the 
invariant density on the attractor, are used as constraints on the choice of mapping 
parameters. We give a detailed analysis of methods to extract the Lyapunov exponents 
and show how to equate them to the values for the parametric map in the constraint 

optimization. 


* If the data points are uniformly distributed in embedding space, our data are truly 
stochastic. But if the data points are distributed in a bounded region (attractor), then 
the data are not stochastic, and order can be extracted from data. 
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If the evolution of our dynamical system is on such an attractor, then the d-dimensional embed- 
ding space enclosing the attractor should be sufficiently larger than d A that all the geometric 
information about the attractor will be exposed in the embedding space. In the next section we 
present the methods of forming this phase-space reconstruction. 

2.1 Takens-Packard Phase-Space Reconstruction 

Using this method, we construct from the solar flux time series x(t)’s d-dimensional vectors 
which, when embedded in R d , describe the full dynamical evolution of the system. Section 3.1 
is devoted to the techniques used to identify the correct value of d directly from the time senes. 

Suppose we have the correct value of d. We consider measuring a single scalar variable x at 
discrete time points x(n) for n= 1,2, . . . ,N 0 . In our case, x is the solar flux and the observa- 
tions are daily. 

To illustrate the technique, we have exaggerated the sampling intervals and magnified a small 
portion of our time series in Figures 2-2 and 2-3. Here q and q should be extremely close to 
one another and r should be chosen so that no information is lost. The choice of r will be dis- 
cussed later. It cannot be too big or too small. For right now, suppose we have the correct value 

of r. 

Figure 2-2 shows hov^can construct vector kets |y x >, |y 2 ), ••• IYn)- 

The set of | y(n)>, of which we have N = N D - d, captures the evolution of our nonlinear dynam- 
ical system as it moves through the d-dimensional phase space. Familiar phase-space coordinates 

are time derivatives x(n), x(n), x(n) evaluated at discrete times. The times, lagged x(n), are 

nonlinear combinations of the local time derivatives and are fully acceptable substitutes for the 
usual phase-space coordinates, as Eckmann and Ruelle have emphasized (Reference 7). 
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Having the y(n) and the embedding space, we ask how we can use the series of yin) to predict 
|y(N+l)), | y(N+2)>, etc. That is, given a data set |y(l)>, | y(2)>, ... |y(N)), is there a mappmg 


F from R d to itself parameterized by | a> = 


al 

al 


ap 


that takes us from |y(n)) to |y(n+l))? 


|y(n+l)> = F(|y(n), |a>) 

that is |y(2)> = F(|y(l)> , |a>) and |y<3)> - F(|y<2)>, |a» 


For simplicity we write F = F and | y > = % I a > = S and generate tbe columns in Figure 2-4. 

Therefore, from a time series we can form phase-space reconstruction in embedding space of 
dimension d, as shown in Figure 2-5. 

Our function F(J,S) comes from parametrically "fitting" the right-hand column in Figure 2-4 of 
J(n+1) resulting from the left-hand column of ftn). Fitting the data then suggests making an 
estimation of 2 so that a cost function C(2) is minimized. 

/ 

N - 1 d (2-2) 

C(S) - Y, E (y>-i)-F,[K")^) 

/i-l \m m 1 - / 

Here we are not just making a fit to data with a set of functions F(ft2). Rather, these functions 
evaluated along the orbit are to be related to each other in the manner of a dynamical system. 
That means the function F(?,2) evaluated on the data vector %n) is required to do more than 
reproduce y(n+l). F(y,2) must also be a function that, when iterated, will reproduce y(n+2) 
after two applications to ftn), ftn +3) after three, and so on. That is, the geometrical properties 
of our dynamical system given by F(ft2) are used to determine the success of the fit. 

The data contain invariant dynamical information that is essential for a full description of the 
structure of the attractor that they evolve on. We would try to impose these invariants as con- 
straints on the fit. This emphasizes that we are creating a dynamic and not just a fit to data. 
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Figure 2-4. Representation of a Map 
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Figure 2-5. Taken-Packard Transformation 
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Many invariant dynamical quantities remain constant as the system evolves m time. We will 
concentrate on two kinds that allow us to construct a reasonable predictor. One kind of invariant, 
the Lyapunov exponents X*, ... X d> describes the expansion or contraction of phase-space 
volumes under the iteration of F(y,2). The second kind of invariant is the density of points on 
the attractor p(y). It captures global features of the frequency with which orbits visit various 
portions of the attractor. It is a different kind of invariant from Lyapunov exponents. Its in- 
tegrals with smooth functions G 0) are unchanged under operation with the mapping function that 

underlies the dynamics y(n) -* y(n+l). 

fd pCDGCD = fd d yp(y)G(F(y,a)) (2-3) 

Both the Lyapunov exponents and invariant density are independent of the initial conditions on 
the orbit. 

Here we find the parameters IT in F(y,2) minimizing a cost function subject to certain constraints. 
The constraints are chosen so that iterations of the mapping function F(y,2) give nse to values 
of dynamical invariants that are the same as those indicated by the experimentally measured data 
set y(n). In this way essential geometric information about the particular attractor on which the 
data are found will be built into the parametric mapping. In the next section we discuss the 
structure of the predictor. 

2.2 The Structure of the Predictor 

Now that we have successfully embedded the data x(n) in R d by creating d-dimensional vectors 
y(n), n=l, N, we need to choose a class of parameterized mappings, a cost function to 
minimize, and a means to impose the constraints on our minimization. Our maps are required 
to look around at the behavior of the phase-space neighbors of the point ?(n) and predict forward 
according to how a cluster of phase-space neighbors, regardless of its temporal sequence, is 
moved forward in time. The map will then map forward any new point y to some weighted 
average of its neighbors’ forward evolution. 


10000078 


2-9 



We take our mappings to be of the form 


F(y,a) = £y(n+l) g [y,y(n);a ] (2 ' 4) 

n-1 

where g[y,y(n); y] is near 1 for 7 = y(n) and vanishes rapidly for nonzero | y - 7(n) | . This type 
of mapping is very similar to the form used in communications engineering under the name 
"kernel density estimation. " 

There are a few general requirements that g[3, 7(n); 3] has to satisfy. These requirements are 
common in estimation theory (e.g., as cr 0, g[y,7(n);2]-*‘ delta function). The value of 
gfy.?( n );^] should be numerically stable and computationally efficient. These requirements are 

all satisfied by the following choice. 


g(y>y(.n);a) = 


- 1> - j^rt) I 2 

e 

P 

a . + ay(n)(y - y(n)) + £ a k 

i-3 

\y - Kn)| 2 ] 

° J 

m k 

-\y - yt ")! 1 

P 

£-3 

' \y - K»)l 2 ' 

J 

m k 



The parameter 3 is P-dimensional, 3 = [a^a^. . . .a 7 ] + , and a is a fixed parameter for scaling. 
The variable m t is also fixed at various values. We can treat a and m t as parameters to be 
optimized in the same sense as the a. 

The choice of the cost function is also up to us. Since we treat £(7,3) as a dynamical system 
evolving points 7(n) into y(n+ 1), our map should reproduce from 7(n) not only 7(n+ 1) but, via 
application of the map, a sequence of 7(n+l), 7(n+2), . . . y(n+L) U P t0 some L beyond which 
we simply do not trust the accuracy of our algorithm F or of the computer we use to compute 

the future ys. 
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In electrical engineering it is a common practice to form predictors of linear systems as 

y(m + 1) = 52 %k ~ ^ + ^ ^ ^ 


This suggests that the predictor for future points of a nonlinear system is a linear combination 
of iterated powers of the map F(^,3) as a generalization to the linear predictor. The nonlinear 

generalization is 


y(m + 1) - 52 X K F k [ y(m - k + l),a ] 

i-i 


(2-7) 


where F k is the kth application of F. If Z<$& were the exact mapping, we would require 
Y^X k = 1. The X t weight the various iterations of F and are used to determine which iterations 

°f(£jvve believe are the most accurate. We should require Xj > Xj +1 to indicate that the lower 
iterations of F are believed to be more accurate than higher iterations. 

Our predictor combines both past information from times m-k+1 and k = 1,2,...L and infor- 
mation from all the phase-space neighbors of the orbit points y(m-k+l) because of the structure 

of E(y,a). 

By extracting the phase-space information in F(y,2), we efficiently tap properties of the full data 
set. In the next section we study the form of cost function for our predictor. 

2.3 Cost Function for Nonlinear Predictor 

The cost function is nothing but a normalized rms deviation from our prediction to the actual 
data points. Therefore, the cost function associated with our nonlinear predictor is 


C(X,a) = 


N - 1 


n=l 


L 

y(« + i)-E^i pk + M 

k* i 


( 2 - 8 ) 


52 iw ■ yoof 2 

rt-i 
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This kind of cost function will automatically contain information on the Lyapunov exponents, 
which themselves are expressions of the dynamics as iterations of the map. Some information 
on the invariant density function on the attractor is also contained in this unproved cost function. 
Different choices for the function g[y,?(n); 3] can be taken. The Gaussian we work with could 
be replaced by a Lorentzian or other choices that weight neighbors more. 
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Section 3 r Invariants of Dynamical Systems 


Having a map and a cost function, we are ready for the constraints. We discuss how to deter- 
mine Lyapunov exponents and invariant density from the E(J,S). Equating the numerical values 
for the Lyapunov exponents X, extracted from data to their expression found from E(J,D will 
give us our fust set of constraints on the minimization of £<*,S). These constraints produce the 
optimum values of parameters 2. 

Since we have a finite amount of information, we choose to express this in terms of the projection 
of p(50 on a set of dual basis functions that is a complete set in R d . By projecting the p($) 
determined from the data onto these basis functions, we can determine the coefficients of the 
expansion of pQ) in this basis. Similarly, we can project the p(y) determined from the map 
F(y,a) onto these basis functions and determine the expansion coefficient of the map. Equating 
the coefficients from the data to the ones determined from the map gives our final constraints 

on the minimiz ation of C(X,2). 

One can test these techniques by numerically generating a data set of x(n) from the known 
attractors like Henon or others and treating these data as having come from an unknown source. 
As the dimension of the phase space increases, the amount of data necessary for accurate 
prediction increases dramatically. Once F 0,2) is found, the details of F(?,2) should reflect the 
known features of the phenomena giving the signal. 

3.1 Choice of the Embedding Dimension d 

Here we would like to determine the correct value of the embedding dimension d from the scalar 
time series x(n), n=l,2,...,N D . We assume that there are enough data that we need not be 
concerned with statistical issues about numerical accuracy. We also assume that extrinsic noise 
is absent from the x(n) when we receive them. We further assume that by following Taken’ s 
phase-space reconstruction technique we have successfully captured the dynamics and embedded 
our time series. This requires a correct choice of r, which will be discussed in the next section. 
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For now, let’s further assume we have a correct r to construct the attractor in the phase space. 
To establish dimension d, we need some characteristic of the attractor that becomes unchanging 
as d becomes large enough, thus indicating that the attractor can be embedded in R d . This 
invariant characteristic of the attractor is the attractor dimension d A . One increases d until d A 
remains constant and identifies the minimum d where d A "saturates" as the embedding dimension. 
But computation of d A is difficult, so we use the correlation function D(r) proposed by Takens 

(Reference 8). 

dwm> - £ £ U(r ' m ' K0I) ‘ * 1 <M) 

Jl z>0 

where U(z) is just the unit step function U(z) = z<0 • 

For N large enough, the behavior of D(r,N,d) for r becomes independent of N and D(r,N,d) takes 
the form 

D(rfl,d) = *( rAr** ® ^ 

If we plot D(r,N,d) versus r we can single out the correct value of dimension d as in Figure 3 1. 

From Figure 3-1 it is concluded that the minimum value of d=3 is the right choice beyond which 
attractor dimension d A does not change or the slope of our graph becomes constant. 

In the next section we study the correct choice of r to reconstruct the phase-space attractor. 

3.2 Choice of the Time Shift r 

The choice of time shifts r is not well agreed upon. If the underlying system were a differential 
equation and a scalar variable x(t) were measured at discrete times x(n) = x(t 0 + n At), then 

we would be, by the choice of logged variables, trying to find a discrete replacement for the usual 

phase-space coordinates: 

d d ~ l x(t) 

x(t), 


10000078 


3-2 



I 


I 



Embed- 



The best choice for time shift r is a fraction of the time associated with the first local minimum 
of the autocorrelation function 

— f T x(t + x)x(t)dt 

T Jo 

We find that this choice, although somewhat arbitrary, works well in practice and gives a simple 
systematic way of determining r. The autocorrelation of more than 4000 data points of solar 
flux time series is shown in Figure 3-2. The first local minimum occurs at t= 13 days, so r=3 

days is a good choice. 


3-4 
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Aiilocovoilancus 


Autocovariances 



Autocorrelation and its error 



P!ot of Autocorrelation and 
Standard Error of Autocorrela- 
tion for Solar Flux Time Series of 
300 Shifts 


Figure 3-2. Autocovariance and Autocorrelation of Solar Flux Senes 


Plot of Autocovariances for Solar 
Flux Time Series of 1 200 Shifts 
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Section 4. The Most Important Invariant of a Time 
Series 

4.1 Extracting the Largest Lyapunov Exponent from a Time Senes 

t 

4.1 .1 Description Of Extraction 

The sum of the Lyapunov exponents is the time-averaged divergence of the phase space trajec- 
tory; hence any dissipative dynamic system will have at least one negative exponent. Any 
dynamic system without a fixed point will have at least one zero Lyapunov exponent. 

A small positive Lyapunov exponent is an indication of chaos and a very large positive Lyapunov 
exponent is an indication of a totally stochastic or random system. Therefore, the sign of the 
exponent provides a qualitative picture of a system’s dynamics. A positive exponent represents 
chaos, a zero exponent represents a marginally stable system, and a negative exponent, a periodic 

system. 

Figure 4-1 shows the actual solar flux data and their largest Lyapunov exponent for more than 
4000 points. Hem we have used the technique of phase-space reconstruction with delay coor- 
dinates shown in Figure 2-2 for a small portion of a time senes. 

To check whether the program that generates Figure 4-1 is really fimctioning well, we have 
plotted in Figure 4-2 the Lyapunov exponent of a time series that was generated from a sinusoidal 
Omt-tinn Because a sinusoidal function is well behaved and totally deterministic, its Lyapunov 
exponent should approach zero. (Remember that we can generate only a truncated sinusoid and 

not an infinitely extended sinusoid.) 

This technique of extracting the largest Lyapunov exponent was originally developed by Wolf, 
Swift, Swinney, and Vastano (Reference 5). The disadvantage of this method is that it gives only 
the positive largest Lyapunov exponent. To improve the technique that can enable us to 
iattnmc all the Lyapunov exponents, Eckmann, Kamphorst, Rnelle, and Ciliberto developed 
a QR matrix decomposition method of calculating the exponents. Their method is more 
r u^,.H than the Wolf-Swift-Swinney method and is computationally more expensive. For 
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Figure 4 - 2 . The Lyapunov Exponent of a Truncated Sinusoid 
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this reason we refer the reader to Reference 11 for further investigation of their method. 
However, because the Wolf-Swift-Swinney method is very simple, we will use it in our construe- 
tion of a predictor and will discuss this method in ten steps. 

4.1.2 Procedure for Extraction: The Wolf-Swift-Swinney Method 

NPT = length of our solar flux time series 

DIM = dim ension of the phase- space reconstruction 

TAU = time delay reconstruction 

DT = time between the data samples -* to normalize the exponent 

SCALMX = length scale we consider to be larger -*> no more information is probed t > scalmx 

SCALMN = length scale we consider to be small -* noise dominates i < scalmnx 
EVOLV = constant propagation time (A=t k+1 - t*) 

ANGLMX = maximum angular error (not a free parameter): does not really change the 
exponent estimate and is usually between 0.2 and 0.3 radians 

1) Search data for nearest neighbor to the first essential point. (This is not easy.) 

2) Disregard points closer than SCALMN. 

3) Propagate current pair of points EVOLV steps through the attractor and compute the final 
separation. 

4) Find the log of the ratio of final to initial separation of this pair. This updates a running 
average rate of orbital divergence. 

5) Attempt a replacement step. 

6) Calculate the distance of each delay coordinate point to the evolved essential point. 
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7) Check whether or not the angular orientation is less than ANGLMX radians for points 
closer than SCALMX but farther away than SCALMN. 

8) If more than one point is found, use the point defining the smallest angular change for 
replacement. 

9) If no points are found, take replacement points as far as 2 x SCALMX away. 

10) Repeat this process until the essential trajectory gets to the end of the data file, by which 

time we hope to see stationary behavior in the Lyapunov exponent. 


X = 


M 


‘o' 


log* 




(4-1) 


wh ere M is the total number of repl acement steps. 

/''TfvTsZr - A 'ft'"®*'*. ~ 

C A^ Figure 4-1 Thows, we haveTgood convergence of the largest Lyapunov exponent for 

parameter values of N=4090 data points, d=3 three-dimensional embedding space, r =3 three 
day time shift, DT= 1 one day time difference between the samples, SCALMX - 0.05 maximum 
length scale, SCALMN =0.00001 min imum length scale, EVOLV=4,5,10,15. As we increase 
the EVOLV, we reach a point where the exponent does not change its convergence behavior. 
The first value of EVOLV that makes a satisfactory convergence is the optimum EVOLV. In 
our case, this value is EVOLV=4 and the Lyapunov exponent is about >w = 0.01. Now that 
we have the Lyapunov exponent directly from our time series, we are ready to calculate it from 
the map. Setting equal X calculated from the map and X calculated from the time series would 
give us a set of constraint equations to calculate parameters 1 to minimize the cost function 


C(X,2). 

4.2 Determining the Lyapunov Exponents from the Map £(?,§) 


Whatever method we use (the Wolf-Swift-Swinney method or the Eckmann-Kamphorst-Ruelle 
method) to determine the Lyapunov exponents from the time series, we must now establish a way 
to express these same quantities in terms of our map £<*,*). This was originally done by 
Shimada and Nagashima (Reference 12). But to easily use the results in our optimizations, we 
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Figure 4-3 Evolution and Replacement Procedures Used to Estimate Lyapunov 
Exponent Directly From a Time Series 
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can use the slightly different approach presented by Abarbanel et al. (Reference 13). 

Lyapunov exponents characterize the way neighboring points, small areas, or small volumes near 
Ore orbit of interest evolve under the mapping. To find them one linearizes the mapping 
y( n+l) = HJ(n),S] around a given orbit J(l), ?(2), J, . • • ?(■>)■ Small deviations from this 

orbit, called 5?(n), evolve as 

5y(n+l)= M (y(n)] 5y(n) 


where 


[M(y)] jt = Fj<y,a) 

°y t - 


(4-3) 


is evaluated along the orbit of interest. The Lyapunov exponents are found from the eigenvalues 
of matrix M k [y(l)] 


AT* 

y(D 

11 

a: 

A 

-I 

M 

K*-D 

M 

yik-2) 

..M 

y(X) 

- 

. 

- 

- 


■ 


■ 




Now, to find the largest Lyapunov exponent, we apply the Matrix M k to an arbitrary vector # . 
Then, forming 




= X. 



\M\W)\\ 

l|vv& _ 


(4-5) 


from the trace of the Matrix M k we can write 



— Ln [triM*) 
t k 


(for large k) (4-6) 
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Remember that M is a function of 3. 


To find the next largest exponent \ z we can use 
— Ln i[tr(M *)] 2 - 

* ~ Jj- 
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The frequency with which orbits y(n) visit regions of the phase space R J defines an invanant 
distribution function p<?), which is formally defined for the mapping ?(n+l) = F[J(n)] as 


p(y )= ta ’f 1>'-F‘[>'(1)]> = (for maps) (5-1) 

N ~~ N jt-i 

In similar fashion, the invariant distribution for a numerical data set y(n), n= 1,2,...,N is given 
by 


p(y) = 


lim J_ 
v- N 


N 

£ 5 J [y-yW] 
k= 1 


(for time series) (5-2) 


Any finite sequence of N points has a finite resolution on the attractor. That resoluuon is 
approximately N‘ lfdA , which is the order of the mean distance of N points on a d A -dimensional 
set To handle this matter of finite resolution, we introduce a complete orthonormal set of 
functions M ) defined on R a , which can serve as a basis set. Truncadng this expansion at some 
finite order p=G provides a finite-resolution representation corresponding to whatever infer- 
mation we have on p(7). We then expand 


G 

p(y) = E B .* 

\i=i 


(5-3) 


The coefficients B„ will be the invariants of the dynamic process that characterizes pQ) within 
a given basis ^,(y). 

Now, the challenge is to extract B, from the data vectora J(n) and from the parameterized map 
E(J,a). w q.v.H..; the B, from the data to those from the map will produce our final constraints 

on the cost function C(X,2). 
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5.1 Coefficients B, Of The Invariant Measure From The Time Series 

A„ y choice of basis *,® would do the job. such as the complete orthogonal (or orthonormal) 

set used in the Fourier Series «J) = <=*'. where a = [m„ m 2 . m 3 , ... mj. But since our 
attractor is bounded, most of the work performed by the Fourier representation of P® will be 
expended in making p® vanish from the attractor. Thus, wha, we need is orthonormal fimcttons 
concentrated on the attractor. Now we construct an optimal choice for *,®. Take the total data 
set ?(n), n = 1,2, .... and divide it into two portions as in Figure 5-1. 



Figure 5-1. 


Breaking the Time Series Into Small Sections for Proper Treatment 


Now take the second portion of length Q and further divide it into G subsections, each L long 
(Q=GL). Each group is a sample of the invariant attractor. Treat each of the G data sets as 
an independent sample of p® and form the invariant distribution for the crth sample (see Figure 

5-2). 



Figure 5-2. One Portion of the Time Series 
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treated as data lor modeling 





p 


,w/ 


L 


with a = 1, 2, . 
the mean density 


.. G. The data point J(k,a) is the kth member of the ath sample. Of course, 
of the G samples is just the total invariant density of the data set of length Q. 



P«()0 


(5-5) 


Now, from G samples pa(y) we form the following phase-space correlation function: 


R(ZW) . J_ V pjz) p a (w) (Every iteration is within a given sample. (5-6) 

G a- 1 

It can be shown that the normalized eigenfunctions of this correlation function are the optimal 
^(7) for expansion of functions localized on the attractor. 

The requirement that ^(?) be an eigenfunction of R<2,*) leads t0 

(5-7) 

fd d ‘ Rlytd *,(z) = 


The are normalized as follows 

f d d * *,( w )* u '(*0 = 


(5-8) 


For an infinite G, the set of eigenfunctions becomes complete. For finite G, R(W,Z) becomes 
a finit e sum of separable kernels. Thus 


*,(50= E 0.(50 

--i 


(5-9) 


Now they are localized near the attractor, because pJSft vanish beyond the attractor. Simple 
substitutions reduce the problem to a finite matrix problem. 
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(5-10) 


d d z 


-- E p.©p«w> 

v a*l 


E c «p«© 


a - 1 


= H E C «Pa(^) 


a- i 


2 A a C p “ = nC wfcere ^ = ^\d d Z9 a (A?^ 


(5-11) 


p-i 


using the nor maliza tion condition we have 


G / 

E C <X 

a = 1 


-L 5 p-jA 7 

\iG 


(5-12) 


This shows that all the eigenvalues /i are positive. 

For all the calculations and computations of pJS) we use the Gaussian representation of Afunction 

-/ s (|x|) (5-13) 

'J(nu) d 


therefore 


p«(y) 


i L i L 

= T E h d> - **.«))! -T E 


-!?-**,«) I 2 5 


K = 1 


L tl 


(5-14) 


jfc=l 


thus 

= 1 1 y' e -|**,«)-y(/.MI 2 /* (5-15) 

“ p Votw) 5 gl2 «“ l 

when » is small, we have only a small loss of resolution in calculating pj$). Now we are ready 
to find our optimal fJS) from the G data sets. First we calculate the G x G matrix A*. Next 
we calculate the eigenvalues p and eigenvectors CJ of this matrix, being sure to normalize them 
according to Equation (5-12). Now we can calculate i/ffl by using the normalized Q and the 
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ptf) from Equations (5-9 and 5-14). A typical value for L is 750 ? points and for G is 5 
samples. These are good values for N=4500 Xpoints of our time series. We now have G 
orthonormal functions ^(7) extracted from G samples p(7) of the invariant distribution. We can 
use Equation (5-8) to project a particular B„ out of Equation (5-3), 

= f d *y p (y) 4^,00 (5 " 16) 


This shows that B„ are invariants of the dynamics since they are integrals of ^ with the density 

p(7)- 


Therefore, 


bA± t c: 


N: 


■1 K-l 


G L 

EEw, 


-1 


(5-17) 
(from data) 


because from Equation (5-2, 5-9, and 5-16) we have 

$ r =(d J yp(y)<VJy)‘jd d y E 


E Cp.W 


«*1 


G N 


= ff cl -nE E c - P.W*)) 

JV-«> Ai rrri J rv B =i K =i 

P.6W) 


a=»l ic=l 


(5-18) 


Thus, Equation (5-17) is used to calculate B„ from the data, where B, are the G numbers 
characterizing the invariant density p(7) by its projection on the optimum basis vectors US)- 
In the next section we discuss how one can find B, from the map. The G equalities between 
these two evaluations of B, form our final constraints on the minimization of the cost function 

C(X,J). 
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5.2 Coefficients Of The Invariant Measure From The Map 

To determine B p from the map F(?), we look at the definition of invariant density as expressed 
by Eq ua tion (5-1) and Equation (5-2). Define A K the projection on '/'iAJ) 

A K (|x) = (D)) = (KD)) (5-19) 

This is the projection of 5 d (y-F x (y(l))) onto the orthonormal ^„(y). That is, we can expand the 
5 function in terms of ^(J) to get 

i'b-F'Wil- (5 ' 20) 

H-l 

Therefore Equation (5-1) can be written as: 


1 IV K-1 u-1 H-l 


N 


«K-l 


^(50 


(5-21) 
(for large N) 


Comparing this with Equation (5-3), we have 


V b £ i E f^wi))) 

" K -l iY K-1 

From Equation (5-9) and the definition of p($) we can write this as 

l=if f fc; «<(«/, a) - PW 1») 


L y-i «-i 


( 5 - 22 ) 


(5-23) 
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We can replace F(y(l)) by F(yOt)) and therefore 


i N 

B =—T 
d NL ^ 


K*1 




1 N 
— E 

NLti 


a*l 7*1 




ee ^ 

a-l j-l \J(TZ(j>) d 


(5-24) 
(from a map) 


We also observe that 


b d (y-F K * 


Ky(l))) = 


d d wd d (y-F(w))6 d (w-F K (y(l))) 


(5-25) 


thus 


where 


and 


a k+ 1 (p) = E V' 



c 1 L 

d d yp L (jO^(y)= P *00 = 1 E 4 c(P) 


which is 


(5-26) 


(5-27) 


(5-28) 


(5-29) 
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since 


i> 


lim p L (|i) = 


^(y) = Bp 


(5-30) 


This shows that the B,, are the components of the eigenvector of T with eigenvalue unity, if 
Equation (5-1) for p($) converges. 
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Section 6. Conclusion 

In this document we developed a set of procedures tor forecasting solar tiux without rererence 
to any underlying solar physics. These procedures are applicable to any time series with a 
broadband power spectrum. We introduced techniques of extracting dynamical invariant directly 
from solar flux time series. The existence of positive, small Lyapunov exponents supports our 
hypothesis that solar dynamics is not a stochastic process; it is indeed chaotic. 

A FORTRAN computer program for procedures outlined in this document is under development, 
and forecasting results will appear in future documents or technical papers. 

When we have perfected these programs, we will be able to apply these techniques to any time- 
series predictions (orbit, attitude, atmospheric density, geomagnetic index, ... stock market). 
Perhaps the most important application will be to a time series generated by subtracting the 
predicted positon of a spacecraft using GTDS (with density flag = off) from the observed 
position of a spacecraft from telemetry or ground track data. This allows us to model all 
perturbations on the trajectory of the spacecraft directly from data. 
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